Melting of Lennard-Jones rare gas clusters doped with a single impurity atom 



O 



Nicolas Quesada 
Instituto de Fisica, Universidad de Antioquia, AA 1226 Medellin, Colombia 

Gloria E. MoyancQ 
Instituto de Quimica, Universidad de Antioquia, AA 1226 Medellin, Colombia 

(Dated: July 13, 2010) 

Single impurity effect on the melting process of magic number Lennard-Jones, rare gas, clusters 
of up to 309 atoms is studied on the basis of Parallel Tempering Monte Carlo simulations in the 
canonical ensemble. A decrease on the melting temperature range is prevalent, although such effect 
is dependent on the size of the impurity atom relative to the cluster size. Additionally, the difference 
between the atomic sizes of the impurity and the main component of the cluster should be considered. 
We demonstrate that solid-solid transitions due to migrations of the impurity become apparent and 
are clearly differentiated from the melting up to cluster sizes of 147 atoms. 

PACS numbers: 36.40.Ei, 61.46.+W 
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I. INTRODUCTION 



Alloying effects in atomic nanoclusters cover a domain 
of property behavior wider and more complex than those 
corresponding to individual atoms and bulk matter, with 
strong particle size specificities which combine with com- 
position and finite-size effects. Even for pure substances 
the structure of their atomic nanoclusters is very de- 
pendent upon the number of atoms per particle. There 
are "magic" numbers, corresponding to cluster structures 
characterized by their conspicuous energetic stabilities 
relative to size, but for a given finite cluster structure, 
stability results from a trade-off between packing and sur- 
face effects. General non-monotonic property trends as 
a function of size characterize finite clusters, so complex 
structural transitions may occur during the growth from 
finite sizes to the bulk. The addition of dopant atoms to 
a pure atomic cluster can alter its structure and growth 
patterns depending upon the nature of both the impurity 
and the cluster, the cluster size, and the concentration of 
dopant atoms. The possibility to manipulate nanoparti- 
cle structures and so, tune their physico-chemical proper- 
ties {e.g. catalytic, electronic, thermodynamic) has mo- 
tivated a lot of recent research on alloy nanoclusters^,. 
Regarding the phase changes, the melting process of pure 
and alloy clusters has attracted considerable attention in 
experimental as well as in theoretical studies. A number 
of specific features have been recognized in the melting 
mechanisms of finite particles such as solid-solid structure 
changes prior to melting 3 -, premelting^ effects of surface 
loosening (formation of "liquid- like" surface layers 4 ^), 
coexistence of different atom-packing schemes^, oscilla- 
tions between the liquid and solid phases^, etc. 
The melting temperature as a function of the cluster 
size has been studied on the basis of several models 
which agree on predicting that the melting temperature 
decreases linearly or quasi-linearly with the inverse of 
the radius of the particle^^. The pioneering work by 



Pawlow is summarized in the formula: 

T M (A0 = T M (co)(l - C7V- 1 / 3 ), 

in which Tm(N) and Tm(oo) represent the melting tem- 
peratures of a A-sized spherical cluster and the bulk, 
respectively; and C is a constant (see Ref Q for a deriva- 
tion of this law and further correction terms). Pawlow's 
law is consistent with several experimental results and, 
although deviations occur for the smaller clusters, whose 
shapes are far from spherical, the melting point of nan- 
oclusters is usually depressed. Nevertheless, there is ex- 
perimental evidence of exceptions to this trend for cases 
like the ionic tin clusters with 10-30 atoms, whose melt- 
ing points are at least 50K above that of the bulki£. In 
addition to the size effects, the melting temperatures of 
alloy clusters can be increased 1 ^ or decreased 1 ^ with re- 
spect to those of the pure components 1 -. The amount 
and direction of the shiftings of the melting point in fi- 
nite doped atomic clusters can be attributed to several 
factors: alterations of the cluster structure, whether or 
not the impurity is soluble in the cluster, many-body en- 
ergetic effects, and/or other complex energetic-entropic 
effects 1 ^. 

The phenomenology seen in the melting mechanisms of 
pure clusters is also apparent for binary and multiple- 
component clusters, but having the composition as an 
additional variable enormously increases the complexity 
of structural behavior a 13 : 14 . Alloying effects in mixed 
atomic clusters depend upon the differences between the 
atomic sizes, cluster surface energies, overall structure 
strain, number and strength of the interactions between 
unlike atoms^ 4 . Further contributing aspects may be 
kinetic factors, specific electronic/magnetic effects, and 
environmental conditions 1 . Alloying effects can be sig- 
nificant even when a single impurity is introduced into a 
cluster of the order of a hundred atoms^. 
An efficient scheme to model the melting of doped atomic 
clusters has to address the issues associated with the in- 
creased complexity of the energy landscapes to explore 
during the simulations of mixed clusters, the occurrence 



of homotop structures, as well as convergence difficul- 
ties related to quasi-ergodicity that have been described 
elsewher o 15 ' 16 . Methods such as replica exchange Molec- 
ular Dynamics and Parallel Tempering Monte Carlo 
(PTMC) have been developed to address the quasi- 
ergodicity by improved sampling. PTMC is a powerful 
method to sample rugged energy surfaces which takes ad- 
vantage of the fact that replicas running at high temper- 
ature are able to sample most of the relevant configura- 
tion space. At the same time, through configuration ex- 
change PTMC connects high temperature replicas, which 
can visit most of the configuration space, with replicas at 
low temperatures so that the latter do not get trapped 
in local minimal. 

The paper has been written as follows: In section II we 
present the methodology for optimal structure search, 
sampling and observable calculations to monitor the clus- 
ter melting process. Then, in Section III we discuss the 
features that differentiate the melting of doped clusters 
from that of the pure ones, taking into account their com- 
position and cluster size. Special detail is given to the 
study of the low temperature solid-solid transitions. Fi- 
nally we present some general conclusions. 



II. METHODOLOGY 

In this work we used the scaled Lennard- Jones (LJ) 
parameters <n and ej for the rare gas interactions reported 
in H. 
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A. Optimal Structures 



TABLE I: Global minima for the different compositions con- 
sidered, corresponding energies in absolut units E°/eAr, in 
units relative to their matrix composition E°/ei, shell posi- 
tion of the dopant atom in the structure (number of shell 
containing the impurity/total number of shells in cluster, the 
th shell is the geometric center of the icosahedron) , point 
group, and the melting temperature of each cluster in Kelvin. 



To obtain the (putative) global minima presented in 
Table U (excepting for the cases of the pure L J clusters, 
the 13 atom clusters, Ar^Xe and ArXe54 which had al- 
ready been reported in [14]fl8| ) we perfomed three types 
of calculations: 

1 . Local optimizations using the Fletcher- Reeves con- 
jugate gradient algorithm (FRCGA) were per- 
formed starting from the structures of the global 
minima of each pure cluster, in which one atom 
of the pure cluster was substituted by the dopant 
atom. This way we obtained a set of icosahedral 
low energy structures. 

2. In a complementary, ampler search, we used the 
Basin-Hopping method(BH) 19 . To sample the en- 
ergy surfaces two types of random moves were per- 
fomed: Moving all the atoms at the same time 
and swapping the dopant atom with an atom of 
the matrix. We performed at least 20000 steps 
(=swaps+moves) in which, after each move, we 
performed a local optimization using the FRCGA. 
For all the compositions the BH method arrived to 
the same result of the first procedure. 



3. Additionally, after the finite temperature simula- 
tions described in section Hi B I we quenched samples 
saved at different temperatures for each composi- 
tion. For the smallest clusters we performed around 
25000 local minimizations, and for larger clusters 
about 55000 local optimizations. 

The results were equivalent for all procedures in the 
above list. We note that the first strategy was computa- 
tionally much cheaper than the other two. The minima 
in Table U were used to initialize the finite temperature 
simulations. 

Considering the lowest energy structures the dopant 
atom takes the central position of the cluster when the 
impurity is Ar or Kr, while it remains in one of the two 
most external shells when the impurity is Xe. 



B. Sampling Strategy 

To sample the complex energy surfaces of our systems 
in the Canonical Ensemble we used the PTMC method 15 . 
For each replica we have used two types of moves. On the 
one hand, single particle moves (SPM) have been imple- 
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TABLE II: Cluster sizes (N), number of temperatures sim- 
ulated (n), their minimum (To) and maximum (Tf) values, 
constraining radii (R c ), number of Monte Carlo steps (Nmc) 
and frequencies at which swaps between adjacent replicas were 
attempted (N 3map )- The constraining radii R c and the tem- 
peratures To and Tf are given in units of the LJ parameters 
of the atoms of the matrix. 



merited using an adapstive step that assures that half of 
the time the new configuration will be accepted. On the 
other hand, since we have two different atomic species in 
each cluster we have also implemented particle exchange 
moves. This sampling strategy consists in exchanging 
the position of two different atoms in the clusters. The 
simulation temperatures were chosen according to the ge- 
ometric progression Ti = TqX 1 . The number of tempera- 
tures for the simulations as well as their maximum and 
minimum values are summarized in Table [TTJ 

For each system the number of equilibration steps was 
always equal to the number of Monte Carlo steps (Nmc)- 
To prevent the evaporation of the clusters we imple- 
mented hard sphere constraining potentials for the con- 
straining radii listed in Table [TTJ 

Finally, the swapping acceptance ratios between replicas 
in all the systems simulated remained around 60-70% and 
never went below 35%. 



C. Observables 

We analyse the melting process by monitoring various 
observables. Firstly, the heat capacity Cy, which is cal- 
culated according to the formula: 

Cv(T) = I ^({E%-(E) 2 T ). 

To interpolate the points obtained with the PTMC simu- 
lation and have a smooth dependence in the Cy (T) curve 
we used the multihistogram metho d 20 ' 21 . Note that the 
formula given above depends on the volume in which the 
system is constrained to move. In figure [1] we compare 
the Cy(T) curves for two constraining volumes. No- 
tice that although the second volume is twice the first 
(V(R C = 4.5<JAr)/V(R c = 3.5a Ar ) ~ 2.13) the main 
peak is not strongly affected and the features of the curve 
below the main peak basically do not change (As one ex- 
pects from a "solid" phase). 

To monitor the effects of the dopant atom in each 
cluster we calculated the radial distributions functions 
(RDF) g(r) of the dopant atom and of the rest of the 



atoms in the matrix, for these calculations all the dis- 
tances r have been taken with respect to the geometric 



center of the cluster 



where 
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represents the total number of atoms in the cluster. 
To further quantify the derealization of the atom we cal- 
culated the standard deviation of the RDF of the dopant 
atom (£) according to: 
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D. Harmonic Superposition Method 

To understand the solid-solid transitions that occur in 
a doped cluster between homotops of the same stoichiom- 
etry we have used the Harmonic Superposition Method 
(HSM^. 

This method assumes that there is a number m of 
well defined states that make most of the contribution 
to the partition function in a certain range of tempera- 
tures. Then, one approximates the contribution of each 
state to the partition function (Z(T)) as the contribution 
of its harmonic part. Such partition function is obtained 
from the normal modes and frequencies by expanding the 
potential around the corresponding minimum in a power 
series up to quadratic order: 

V(R) = V(R ) + ]-R T HR + 0(R 3 ), 

where R — (fi,... ,tn), Ro is the equilibrium position 
and H is the Hessian Matrix of that minimum. To obtain 
the partition function (and the thermodynamics of the 
system) one adds the Simple Harmonic Oscillator parti- 
tion functions of each state: 
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where /3 — 1/fcsT, E a is the energy of each state, n a is 
its degeneracy due to symmetry (n a = 2p\(N—p)l/h a , N 
the number of atoms, p the number of impurities, p = 1 
and h a the order of the point group of the state a). v a is 
the geometric mean vibrational frequency of each state 
(which is proportional to geometric mean of the square 
roots of the eigenvalues of the matrix H) and N is the 
number of atoms considered. 



III. RESULTS AND DISCUSSION 

A. Size dependence of the melting temperature 

In figure [2] we present the results of our calculations 
regarding the variation of the melting temperatures as a 
function of N~$ . It is seen that as the size of the cluster 
is increased the melting temperature of the clusters also 
increases, this behavior has been verified for Ar and Xe 



clusters^. The case of N — 13 is certainly out of any 
linear tendency for all the compositions, yet, for other 
larger clusters N = 55, 147, 309 where surface effects are 
less marked the dependence of the melting temperature 
as a function of N 1 / 3 can be well described by a line, 
and in all cases increases with the size of the cluster. 
From figure [2] is clearly seen that doping effects are very 
strong for small doped clusters (N — 13,55), whose 
atoms have the highest differences between their LJ 
parameters, e and a (In this case Argon and Xenon). 
It is also seen that for the largest cluster sizes studied 
here their melting temperatures are almost equal for the 
doped and pure clusters. 



B. Doping effects: Composition and Size 
Dependence 

So far we have discussed the doping effects solely in 
terms of the position, on the temperature scale, of the 
peak associated with the melting of the cluster. Yet as 
is seen in figures [3] and 0] the peak changes, not only its 
position but also its height and width, for some compo- 
sitions. For instance, for Ar^Xe, the change in height 
with respect to Ari4 7 is noticeable although the displace- 
ment of the maximum is just around 1%. Other char- 
acteristic of the Cv(T) that is modified by the pres- 
ence of the dopant atom is the occurence of a small 
peak or bump in the low temperature region. As we 
will demonstrate in section IIII C[ for the clusters with 
sizes (N = 13, 55, 147), this is due to a solid-solid transi- 
tion. Some general trends for compositions ArXe^v-i and 
KrXejv-i (N = {13,55,147,309}) are: Regarding their 
lowest energy configurations, each pure Xejv cluster and 
the doped ArXejv_i and KrXeAr_i clusters have the same 
symmetry group (Ih), and are also very close in geome- 
try. After a small temperature increase, the dopant atom 
in both ArXejv— l and KrXeAr_i behaves the same way. 
It starts to move from the the center of the cluster in the 
lowest energy configuration, to the second most energet- 
ically favorable position, in the outer shell of the cluster, 
as seen in the first two rows of figures [5l El [7] and [8] A 
pictorical representation of the process is given in figure 
[TOl Nevertheless, excepting for the smallest cluster size, 
N = 13, the dopant atom never relocalizes completely in 
a stable configuration different from the global minimum, 
this occurs because for larger structures N > 13 there is 
more than one icosahedral stable structure in which the 
dopant atom is located in the outer shell of the clusters. 
To support this, see in figure [9] that, excepting for the 
cases ArXei2 and KrXei2, the standard deviation of the 
position of the dopant atom (£xe) is always an increas- 
ing function of the temperature, until the cluster melts. 
The bottom rows in figures [5j [6J [7] and [8] show that upon 
melting, the RDFs of the matrix and the dopant show 
the same structure. This indicates that, in the liquid 
phase, and for the compositions studied, Ar and Kr are 



not segregated by the Xe matrix. Finally, as one would 
expect based on the similarities of their LJ parameters, 
the Xe-Kr doped clusters show more resemblance to the 
pure cluster in their Cy curves. 

The clusters Krjv-iXe are the ones that show a more sim- 
ilar behavior to the pure clusters LJjv, considering their 
Cy curves. For these compositions the standard devia- 
tion of the position of the dopant atom (£xe) is always 
an increasing function of the temperature (see figure [9]) . 
This implies that the Xe atom does not leave completely 
its external shell location, as in the lowest energy con- 
figuration (see the fourth column on Tabic [IJ). Such con- 
figuration plays a significant role in the thermodynamics 
of the system until the phase change. This can be seen 
on the RDFs of Xe in the Kr^-iXe clusters, as plotted 
in the fourth column of figures El [3 and [5] Also, on 
the spectra of quenched energies of Kr^Xe and K^Xe, 
in figures [TT1 and IT2l The shape of £xe for Krjv_iXe is 
qualitatively different depending on the cluster size N, 
indicating that the temperature ranges for the migration 
of the dopant atom and the melting of the cluster overlap 
for the smaller sizes. For N — 13, £xe simply increases 
once the cluster starts to melt. For N = 55, the dopant 
atom starts to delocalize smoothly between the second 
and first shells, until the migration is met by the melt- 
ing of the cluster (see the last column on figure |6|) . For 
the largest structures Kr^Xe and Ki^ogXe the dopant 
atom migrates to several positions in different shells of 
the structure, as seen in the last column of figures [7] and 
[U Upon melting, these compositions show the same be- 
havior observed in ArXe^r^i and KrXeAr_i , i.e. there is 
no segregation between the Kr atoms and the Xe atom of 
the cluster. The composition that shows more features 
during the heating process is ArAr_iXe. The largest dop- 
ing effect is seen in the cluster Ar^Xe. For this cluster 
we see that the melting temperature (taken as the po- 
sition of the maximum in the Cv(T) curve) drops by 
around AT — 0.037eA r /kB, which is around 13% of the 
melting temperature of the pure cluster. A comparable 
change in the melting point occurs for A^Xe with re- 
spect to Ar55. This is not the only feature that changes 
drastically when replacing one atom, with respect to the 
melting of Arss . From figure [3] it is also seen that the 
melting peak in the Cy{T) curve for the doped cluster 
is smaller, by almost a factor of 2, as compared with the 
pure cluster, in other words the latent heat associated 
with the melting is smaller in the doped cluster. The 
reduction in the latent heat is a feature present in all 
the Argon clusters, doped with Xenon. For the case of 
Ar54Xe two different transitions are seen in the RDF, 
g(r), of the dopant atom (see the third column of figure 
[6]). These transitions are seen in the non- monotonous 
behavior of the standard deviation of the position of the 
Xe atom in figure [5J3. In the first transition the Xe atom 
migrates from the outer shell to the inner shell, and re- 
mains there. As it was mentioned in the last section, this 
causes a small bump in the Cy curve. Then, as the tem- 
perature is further increased, the atom starts to migrate 
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TABLE III: Parameters used for the HSM. A/t Ar is the en- 
ergy difference between the lowest energy structure and the 
second stable structure considered. A/e t is the scaled differ- 
ence of the two levels considered in terms of the e parameter 
of the matrix. The fourth column is the point group of the 
low energy minimum. Vi is the square root of the geometric 
mean of the non zero eigenvalues of the Hessian Matrix H 
for each stable structure and ksT h = ,„.. _■ , — =4-= — , u n. 

" " (3N — 6) m vq/v\ — In ho/hi 

is the temperature at which the transition between the two 
solid structures is expected, i.e. at which the partition func- 
tions associated with each minimum are equa^. 



tive sampling frequencies of the set of minima obtained 
for each composition are presented in figures QT] and [T2] 
From these figures we note that the extra peak correlates 
extremely well with the appearance of a second stable 
structure that becomes increasingly important until the 
cluster melts. This second structure corresponds to an 
icosahedron in which the dopant atom swaps positions 
with an atom in a different shell. 

To further support our conclusion we have performed 
Harmonic Superposition Method (HSM) calculations for 
some compositions. The input values used in the HSM 
calculations are shown in the first six columns of Ta- 
ble IIII1 In the last column of the same table, we show 
the results (i.e. the predicted temperature for the solid- 
solid change asociated with the transition between the 
two minima ). The predicted temperatures agree well 
with those obtained from the PTMC simulations in fig- 
ures [31 and HI The Table IIIII also shows why the extra 
peak is not present in all clusters. For Kr^Xe, Ar^Xe 
and Kr 54 Xe, as can be seen in the insets of figured the 
temperature of the solid-solid transition is so close to the 
melting peak, that when the structure can change to a 
different minima it has started to sample "liquid like" 
configurations. 



between the center of the cluster, the first shell and the 
outer shell. This occurs near the temperature range for 
the phase change. Finally, when the cluster reaches the 
liquid-like phase an interesting effect occurs, namely the 
Xe atom is segregated from the Ar atoms. This is clearly 
seen in the last row of figures [5J HI [7] and [5] We note that, 
for all the cases studied, the segregation is related to a 
maximum size contrast between the impurity and other 
atoms in the cluster. 



C. Low T behavior 

One of the most interesting features of the Cy cal- 
culations presented in figure [3] and 2] is the occurence 
of a second small peak, not seen in the pure clusters, 
for some of the doped structures. The most noticeable 
case being that of KrXei2- Such peak has been associ- 
ated with a solid-solid transition, and studied in detail 
for rare gas clusters of 6^ and 13^ _ — atoms. It has 
been suggested^— that this bump is due to structural 
transitions between isomers of the same composition. 
We reach the same conclusion via an analysis of around 
1000 structures, which we sampled, for each replica and 
each composition in clusters with up to 147 atoms. We 
later quenched those structures. The energies and rela- 



IV. CONCLUSIONS 

PTMC simulations for rare gases (LJ) doped clusters 
with up to 309 atoms showed that a single atom impu- 
rity can cause doping effects such as the depletion of the 
melting range (with respect to the pure cluster), and the 
occurence of a solid-solid transition in the low tempera- 
ture range^I. The shifting of the melting range due to 
the presence of the single atom impurity decreases with 
increasing cluster size. In terms of absolute temperature 
it is noticeable for clusters with less than a 100 atoms, for 
instance for A^Xe it represents 3.4 K. Several criteria 
(i.e. Cy curves, radial distribution functions, spectra of 
quenched energies, and HSM) have been used to support 
that a solid-solid transition peak may arise for doped 
clusters with up to 147 atoms. 
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FIG. 6: Radial distribution function for 55 atom clusters. 
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FIG. 7: Radial distribution function for 147 atom clusters. 
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FIG. 8: Radial distribution function for 309 atom clusters. 
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